knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
source(here('common_fxns.R'))Create a species richness map for threatened marine species.
Combine taxa-level trend maps to create maps of number of threatened species overlapping intensifying or deintensifying stressors, per cell.
Plot both by absolute numbers of species impacted in a location and by proportion of species (impact counts / species richness).
This is all done at the level of cumulative impact category, not individual stressors.
Set up the dataframe of species to include. By taxon, load all species range maps, smash down to taxon-level species richness; then combine these for a total species richness.
For each impact category, sum increasing and decreasing stressor locations for count:
calc(fun = sum, na.rm = TRUE)Pull the taxon-level maps for each impact category, for both intensifying and de-intensifying impacts.
str_cats <- c('fishing', 'ocean', 'land-based', 'climate', 'all')
taxon_tr_dir <- file.path(dir_bd_anx, 'taxon_trend_rasts')
### Set up filename stems for this run (non-priority)
tr_map_stem <- 'taxon_trend_.+_%s_%s%s.tif'
### .+ is taxon, sprintf args are stressor category, incr/decr, level (1-3)
outfile_stem <- here('_output/rasters/intens_maps/intens_%s_%s%s.tif')
### stressor category, incr/decr, level (1-3)
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
for(trend_direction in c('incr', 'decr')) {
# trend_direction <- 'incr'
for(level in 1:3) {
# level <- 1
outfile <- sprintf(outfile_stem, str_cat, trend_direction, level)
if(!file.exists(outfile)) {
message('processing ', basename(outfile))
trend_map_pattern <- sprintf(tr_map_stem, str_cat, trend_direction, level)
trend_map_files <- list.files(taxon_tr_dir,
pattern = trend_map_pattern,
full.names = TRUE)
trend_stack <- stack(trend_map_files)
trend_rast <- raster::calc(trend_stack, fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(trend_rast, outfile, overwrite = TRUE)
} else {
message('File exists: ', outfile, '... skipping!')
}
trend_rast <- raster(outfile)
trend_pct_rast <- trend_rast / spp_rich_rast
if(any(values(trend_pct_rast) > 1, na.rm = TRUE)) {
stop('some pct raster values > 1, fix it')
}
plot(trend_pct_rast,
col = c('grey80', hcl.colors(20, palette = 'viridis')),
axes = FALSE,
main = sprintf('Pct species with %s %s stressors lvl %s', trend_direction, str_cat, level))
}
}
}### Set up filename stems for this run (priority)
tr_map_stem <- 'taxon_priority_trend_.+_%s_%s.tif'
### .+ is taxon, sprintf args are stressor category, incr/decr
outfile_stem <- here('_output/rasters/intens_maps/intens_priority_%s_%s.tif')
### stressor category, incr/decr
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
for(trend_direction in c('incr', 'decr')) {
# trend_direction <- 'incr'
outfile <- sprintf(outfile_stem, trend_direction, str_cat)
if(!file.exists(outfile)) {
trend_map_pattern <- sprintf(tr_map_stem, str_cat, trend_direction)
trend_map_files <- list.files(taxon_tr_dir,
pattern = trend_map_pattern,
full.names = TRUE)
trend_stack <- stack(trend_map_files)
trend_rast <- raster::calc(trend_stack, fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(trend_rast, outfile, overwrite = TRUE)
} else {
message('File exists: ', outfile, '... skipping!')
}
trend_rast <- raster(outfile)
trend_pct_rast <- trend_rast / spp_rich_rast
plot(trend_pct_rast,
col = c('grey80', hcl.colors(20, palette = 'viridis')),
main = sprintf('Pct priority-wt spp with %s %s stressors',
trend_direction, str_cat))
}
}Net intensification is the difference between the number of species under intensification and the number under abatement.
str_cats <- c('fishing', 'ocean', 'land-based', 'climate', 'all')
infile_stem <- here('_output/rasters/intens_maps/intens_%s_%s%s.tif')
### stressor category, incr/decr
for(str_cat in str_cats) {
# str_cat <- str_cats[5]
for(level in 1:3) {
incr <- raster(sprintf(infile_stem, str_cat, 'incr', level))
decr <- raster(sprintf(infile_stem, str_cat, 'decr', level))
diff <- calc(stack(incr, -decr), fun = sum, na.rm = TRUE)
diff_pct <- diff / spp_rich_rast
diff_pct[incr == 0 & decr == 0] <- NA
hcl.pals('diverging')
brks <- seq(-1, 1, .05)
clrs <- hcl.colors(n = length(brks), palette = 'Red-Green', rev = TRUE)
plot(diff_pct,
col = clrs,
breaks = brks,
axes = FALSE,
main = sprintf('Net intens, Pct spp, %s strs level %s', str_cat, level))
}
}Because a zero net intensification can result from low intensification cooccuring with low abatement, or high of each, let’s examine the frequency with which various combinations occur. Just examine level 1 as this is the most likely to have overlaps.
str_cats <- c('fishing', 'ocean', 'land-based', 'climate', 'all')
infile_stem <- here('_output/rasters/intens_maps/intens_%s_%s1.tif')
### stressor category, incr/decr
for(str_cat in str_cats) {
# str_cat <- str_cats[5]
incr <- raster(sprintf(infile_stem, str_cat, 'incr'))
decr <- raster(sprintf(infile_stem, str_cat, 'decr'))
df <- data.frame(incr = values(incr),
decr = values(decr),
nspp = values(spp_rich_rast)) %>%
filter(!is.na(nspp)) %>%
filter(!is.na(incr) & incr != 0 & !is.na(decr) & decr != 0) ### we're interested in co-occurrence
### only 57989 cells where increase and decrease co-occur?
df1 <- sample_frac(df, size = 1)
tmp <- ggplot(df1, aes(x = decr, y = incr)) +
geom_bin2d(bins = 15) +
geom_abline(slope = 1, intercept = 0, color = 'red') +
scale_fill_viridis_c() +
scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50, 100, 200)) +
scale_y_log10(breaks = c(1, 2, 5, 10, 20, 50, 100, 200)) +
labs(title = sprintf('Intens vs abating, %s strs level 1', str_cat),
x = 'nspp abating', y = 'nspp intensifying')
print(tmp)
}# install.packages("remotes")
# remotes::install_github("slu-openGIS/biscale")
library(biscale)
library(cowplot)
ocean_a_df <- data.frame(prop = values(ocean_a_rast),
cell_id = 1:ncell(ocean_a_rast)) %>%
mutate(a = prop * (res(ocean_a_rast)[1])^2 / 1e6)
### skip ocean - only one stressor, no spp can simultaneously be increasing
### and decreasing
str_cats <- c('fishing', 'land-based', 'climate', 'all')
infile_stem <- here('_output/rasters/intens_maps/intens_%s_%s.tif')
outfile_stem <- here('figs/intens_biplot_%s.png')
for(str_cat in str_cats) {
# str_cat <- str_cats[4] ### 4 = all
outfile <- sprintf(outfile_stem, str_cat)
if(!file.exists(outfile)) {
message('Creating bivariate plot for ', str_cat, ': ', basename(outfile))
incr <- raster(sprintf(infile_stem, str_cat, 'incr'))
decr <- raster(sprintf(infile_stem, str_cat, 'decr'))
df <- data.frame(incr = values(incr),
decr = values(decr),
nspp = values(spp_rich_rast),
cell_id = 1:ncell(incr)) %>%
filter(nspp > 0 & !is.na(nspp)) %>%
mutate(incr_pct = incr / nspp,
decr_pct = decr / nspp)
### create classes: three of intensification (incr), and three of abatement (decr).
### types: "quantile" (default), "equal", "fisher", and "jenks".
### apparently use "fisher" instead of "jenks" for larger datasets
data <- bi_class(df, x = incr_pct, y = decr_pct,
style = "fisher", dim = 3) %>%
mutate(bi_class = factor(bi_class),
bi_int = as.integer(bi_class))
### using jenks:
# 1-1 1-2 1-3 2-1 2-2 2-3 3-1 3-2 3-3
# 2068375 47836 25942 1058722 14297 1768 257262 2743 99
### using fisher:
# 1-1 1-2 1-3 2-1 2-2 2-3 3-1 3-2 3-3
# 2140894 39477 26031 1005903 12893 1254 248652 1841 99
### Quantiles fails b/c non-unique values (in the function).
# quantile(df$incr_pct, seq(1/6, 1, 1/6))
# quantile(df$decr_pct, seq(.05, 1, .05)) ### little of the ocean is decreasing!
# incr_cutoff <- quantile(df$incr_pct[df$incr_pct != 0], .25)
# decr_cutoff <- quantile(df$decr_pct[df$decr_pct != 0], .25)
# data <- df %>%
# mutate(x = case_when(incr_pct == 0 ~ '1',
# incr_pct <= incr_cutoff ~ '2',
# incr_pct <= 1 ~ '3',
# TRUE ~ 'error'),
# y = case_when(decr_pct == 0 ~ '1',
# decr_pct <= decr_cutoff ~ '2',
# decr_pct <= 1 ~ '3',
# TRUE ~ 'error')) %>%
# mutate(bi_class = paste(x, y, sep = '-'),
# bi_class = factor(bi_class),
# bi_int = as.integer(bi_class))
### table(data$bi_class)
bi_pcts <- data %>%
left_join(ocean_a_df, by = 'cell_id') %>%
mutate(a_tot = sum(a)) %>%
group_by(bi_class) %>%
summarize(a_level = sum(a),
pct_level = round(a_level / first(a_tot), 4)) %>%
separate(bi_class, into = c("x", "y"), sep = "-") %>%
mutate(x = as.integer(x), y = as.integer(y))
bi_class_levels <- levels(data$bi_class)
bi_map <- map_to_rast(data, 'bi_int') %>%
rasterToPoints() %>%
as.data.frame() %>%
mutate(bi_level = bi_class_levels[layer])
# create map
intens_map <- ggplot(data = bi_map, aes(x, y, fill = bi_level)) +
geom_raster(show.legend = FALSE) +
bi_scale_fill(pal = "DkViolet", dim = 3) +
ggtheme_map() +
labs(title = sprintf('Intensification and abatement: %s', str_cat))
intens_legend <- bi_legend(pal = "DkViolet",
dim = 3,
xlab = "+ % intens",
ylab = "+ % abate",
size = 8) +
geom_text(data = bi_pcts, aes(x, y, label = pct_level),
size = 3, color = 'white') +
scale_x_discrete(expand = c(0, 0)) + scale_y_discrete(expand = c(0, 0)) +
theme(panel.background = element_blank(),
plot.background = element_blank())
# combine map with legend
final_plot <- ggdraw() +
draw_plot(intens_map, 0, 0, .8, 1) +
draw_plot(intens_legend, 0.6, .6, 0.4, 0.4)
ggsave(plot = final_plot, filename = outfile, height = 4, width = 6, dpi = 150)
}
knitr::include_graphics(outfile)
}